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ABSTRACT 

The temperature in the optically thick interior of protoplanetary discs is essential for 
the interpretation of millimeter observations of the discs, for the vertical structure 
of the discs, for models of the disc evolution and the planet formation, and for the 
chemistry in the discs. Since large icy grains have a large albedo even in the infrared, 
the effect of scattering of the diffuse radiation in the discs on the interior tempera- 
ture should be examined. We have performed a series of numerical radiation transfer 
simulations including isotropic scattering by grains with various typical sizes for the 
diffuse radiation as well as for the incident stellar radiation. We also have developed 
an analytic model including isotropic scattering to understand the physics concealed 
in the numerical results. With the analytic model, we have shown that the standard 
two-layer approach is valid only for grey opacity (i.e. grain size > 10 pm) even with- 
out scattering. A three-layer interpretation is required for grain size < 10 pm. When 
the grain size is 0.1-10 pm, the numerical simulations show that isotropic scattering 
reduces the temperature of the disc interior. This reduction is nicely explained by the 
analytic three-layer model as a result of the energy loss by scatterings of the incident 
stellar radiation and of the warm diffuse radiation in the disc atmosphere. For grain 
size > 10 pm (i.e. grey scattering), the numerical simulations show that isotropic 
scattering does not affect the interior temperature. This is nicely explained by the an- 
alytic two-layer model; the energy loss by scattering in the disc atmosphere is exactly 
offset by the "green-house effect" due to scattering of the cold diffuse radiation in the 
interior. 

Key words: dust, extinction — methods: analytical — methods: numerical — plan- 
etary systems: protoplanetary discs — radiative transfer — scattering 



1 INTRODUCTION 

Protoplanetary discs are planet formation sites. We observe 
the electro-magnetic radiation from the discs to understand 
their physical conditions, and then, to know the planet for- 
mation, especially, its beginning. Property of the radiation 
from a disc is essentially determined by the structure of the 
disc. The structure follows the response of the disc to the 
radiation from the central star. To interpret the radiation 
from protoplanetary discs, therefore, we should have a ro- 
bust link between the radiation from the central star and 
the disc structure, in particular, the temperature structure. 

The temperature structure of the discs is also essen- 
tial for chemical reactions in the discs. The evolution of 
various molecules in the discs and the exchange of these 
molecules between the gas phase and the solid phase on 
grains will be discussed in detail with the ALMA in near 
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future (e.g.. iNomura et alj|2008h . The condensation front of 
icy molecules, so-called 'snow line', is also determined by the 
temperature structur e of the discs (e.g., ISasselov fc Lecan 
l200d : IOka et al.ll2009h . The location of the snow line is very 
important because it significantly enhances the amount of 
solid materials to make planetary cores and affects the sup- 
ply of water to rocky planets. The temperature structure 
also affects the evolution of the discs themselves. The ac- 
cretion activity in the discs is supp osed to be driven by the 
magnetorotational instability (e.g., ISano et alJl200Ch . This 
requires a certain degree of the ionization of disc materials 
which depends on the structure of the discs. 

A milestone in the research of the vertical temper- 
ature structure of p rotop lanetary discs is the work by 
IChiang fc Goldreichl i| 19971 ) (hereafter CG97). They pro- 
posed a two-layer model consisting of the super-heated 
layer directly exposed by the stellar radiation and the in- 
terior warmed by the super-heated layer. This model ex- 
plained the shape of the spectral energy distribution of the 
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discs very well. The computational cheapness of the ana- 
lytic approach makes the CG97 model very useful to com- 
pare with a large sample of the discs observed. Therefore, 
some att empts to refine the simple model of CG97 were per- 
form e d llChiang et alj 200 it iDullemond. Dominik. fc Nattal 
200 ll: IDullemond fc Nattall2003al : iRafikov fc De Corlell2006l : 



Garaud fc Lin 



2007). 



Despite a great success of the CG97 model, numerical 
simulations of the radiation transfer in the discs showed a 
significant decrement of the equatorial temperature relative 
to the prediction by the CG97 model l|Dullemond fc Nattal 
l2003al ). The internal energy loss by the radiation at a 
long wavelength where the optical depth of the disc is 
relatively small is suggested as a cause of the de crement 
ijDullemond et all 12002 ; IDullemond fc Nattal [2003ah . In the 
CG97 model and refined ones, the wavelength dependence of 
dus t opacity was taken into account in terms of mean opac- 
ity. IDullemond et all (|2002h argued the importance of using 
full wavelength dependent opacity. However, we propose an 
alternative approach in this paper: a three-layer model with 
mean opacity which reproduces the temperature reduction 
quite well. In addition, we show that the two-layer approxi- 
mation in the CG97 model is valid only when the dust opac- 
ity is "grey" which is expected if the size of dust grains is 
larger than about 10 (xm. 

Effect of scattering on the vertical temperature struc- 
ture was not considered in the literature very much. Scat- 
tering of the s tellar radiati o n wa s taken into account 
analytically by Calvet et al.l (Il99ll ) and numerically by 
IDullemond fc Nattal l|2003bl ) who showed that the tempera- 
ture of the disc interior is slightly reduced by the scattering. 
How about scattering of the diffuse disc radiation? As the 
size of grains increases, the scattering albedo increases. In 
particular, the albedo of large icy grains is close to unity 
even for the infrared wavelength. This may affect the disc 
structure significantly. Nevertheless, it has not been exam- 
ined so far. 

This paper discusses the effect of scattering of the dif- 
fuse radiation as well as that of the stellar radiation on 
the vertical temperature structure of protoplanetary discs. 
Since the albedo depends on the grain size, we examine the 
scattering effect as a function of the grain size. Although 
there are 2-D/3-D numerical radiation transfer codes avail- 
able publicly, most of them have a serious difficulty in solv- 
ing the radiation equilibrium in very high optical depth 
(r ~ 10 6 ) found in protoplanetary discs (|Pascucci et ail 
12004 ISteinacker. Bacmann. fc Henning 1 120061) . The RAD- 
ICAL developed by IDullemond fc Turollah feoOO) can solve 
such a problem without any difficulty thanks to a variable 
Eddington tensor method. However, it treats scattering of 
only the stellar radiation. We present a variable Edding- 
ton factor code with both scatterings of the stellar radiation 
and of the diffuse radiation but in a 1-D geometry. We also 
present an analytic model to interpret the numerical results. 
This simple model would be very useful to understand the 
physics determining the temperature structure of protoplan- 
etary discs. 

The rest of this paper consists of three sections; in sec- 
tion 2, we develop a numerical radiation transfer code taking 
into account both of the scatterings but only for isotropic 
case and show the obtained numerical solutions. In section 3, 
we construct an analytic model to interpret the numerical 



solutions and discuss the physical mechanism determining 
the temperature structure in protoplanetary discs. In the 
final section, we summarise our findings. 



2 NUMERICAL RADIATION TRANSFER 
WITH SCATTERING 

Our method is an extensi on of the variable Eddin gton fac- 
tor method developed by IDullemond et al" 

I i|2002l ); we in- 
clude isotropic scattering of the diffuse radiation as well as 
that of the stellar radiation. A disc is divided into many an- 
nuli in which the transfer of the diffuse radiation is treated 
one-dimensionally along the normal axis of each annulus 
with neglecting the radiation energy transport among an- 
nuli. This approximation, so-called 1+1D approximation, 
would be reasonable in an optically thick disc, but not at 
the near of the disc inner edge nor in a self-sh adowing re- 
gion (e.g.. IDullemond. Dominik. fc Nattall200ll ). The radia- 
tion from the central star is separated from the diffuse radia- 
tion and is treated with the so-called grazing angle recipe. In 
this paper, we only consider some single annulus cases in or- 
der to feature the effect of the scattering on the temperature 
structure along the normal axis of the annulus. Therefore, 
we assume a grazing angle a = 0.05 radian throughout of 
the paper. The density structure along the normal axis of 
the annulus is solved to be consistent with the obtained tem- 
perature structure assuming the hydrostatic equilibrium. In 
Appendix A, we describe how to obtain the numerical solu- 
tion of the diffuse radiation transfer with isotropic scattering 
in each annulus in detail. 



2.1 A simple dust model 

We adopt a very simple dust model in order to feature the 
effect of scattering on the temperature structure. The ab- 
sorption and scattering cross sections per unit gas mass at 
the wavelength A are assumed to be 



abs 



and 
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respectively. The single scattering albedo at the wavelength 
A is 



sea 



(3) 



The critical wavelength A c may be related to a typical grain 
radius a as A c = 2na. If we consider a spherical grain com- 
posed of uniform material, the absorption cross section is 
expressed as K^ bs = (3VQf 3B )/(<ip d a), where <3* s is the ab- 
sorption cross section normalised by the geometrical cross 
section ira 2 , pd is the grain material density, and T> is the 
dust-to-gas mass ratio. With the values of D — 10~ 2 (Solar 
system nebula), p<j = 3 g cm -3 (silicate), and — > 1 
(A — > 0), we obtain the absorption cross section for small 
wavelengths as 

0.1 nm^ 
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Figure 1. Planck mean properties of a simple dust model as- 
sumed in this paper: (a) absorption cross section per unit gas 
mass and (b) single scattering albedo for the case with u)q = 0.9. 
We show seven cases of grain size from 0.01 um to 1 cm. 

The scattering cross section can be given by the single scat- 
tering albedo for small wavelengths: loo = «;o ca /( K o bs + Ko ca )- 
In this paper, we consider three cases of Wo = (no scat- 
tering), 0.9, or 0.99. The values of ujo for the last two cases 
may be extreme but such a large albedo is expected for icy 
grains in some wavelengths. 

Figure 1 shows Planck means of the absorption cross 
section and the scattering albedo assumed in this paper as 
a function of the temperature input into the Planck func- 
tion. In the panels, we show seven cases of grain size from 
0.01 |i.m to 1 cm. We note that the absorption cross sec- 
tion and the scattering albedo become independent of the 
temperature, i.e. "grey", when the temperature exceeds a 
critical one which depends on the grain size, corresponds 
to the critical wavelength A c , and is roughly expressed as 
T c ~ 10 3 (1 nm/o) K. In this paper, we do not consider the 
size distribution of the dust grains. Thus, the "grain size" 
of this paper means a typical grain size averaged over a size 
distribution function with a weight. 



2.2 Numerical results: Temperature structure 

We here show the results of the annulus with the radius 
of 1 AU obtained from our numerical radiation transfer in 
Figures 2 and 3. The results with other radii have been con- 



firmed to be the same qualitatively. The gas column density 
is assumed to be 10 3 (i?/AU) _1 g cm -2 , where R is the radial 
distance from the central star. The properties of the cen- 
tral star assumed are the effective temperature T* = 3, 000 
K, the radius R, = 2.0 Rq, and the mass M, = 0.5 M Q . 
Other assumed parameters are as follows: the grazing an- 
gle a — 0.05, the visible fraction of the stellar photosphere 
at the annuli / v i s = 0.5, and the mean molecular weight 
flm = 7/3. 

Figure 2 shows the vertical temperature structures of 
annuli with 1 AU radius with various grain sizes. We take a 
coordinate of the Planck mean extinction optical depth with 
the stellar effective temperature as the horizontal axis. Note 
that the maximum optical depth in each curve occurs the 
equatorial plane. The grain sizes assumed are shown in each 
panel. The solid, dotted, and dashed curves are the cases 
of no scattering (i.e. luq — 0), ujo = 0.9, and ujo = 0.99, 
respectively. 

For no scattering cases (solid curves), the so-called two- 
layer structure proposed by CG97 is confirmed. The dust 
temperature near the surface is enhanced due to the direct 
stellar radiation: "super-heated layer" . The thickness of the 
super-heated layer is well expressed by the Planck mean ex- 
tinction optical depth as Tp x * ~ a — 0.05 (grazing angle). 
The temperature rapidly decreases if Tp x * > a. If the interior 
is optically thick against its own radiation, then, the interior 
reaches the thermal equilibrium and becomes isothermal. As 
the grain size becomes larger, the temperature of the super- 
heated layer becomes lower. In contrast, the temperature 
of the interior becomes higher. The physical reason of this 
phenomenon will be discussed in section 3 with two analytic 
models: the standard two-layer model like CG97 and a newly 
developed three-layer model. Here, we just mention the fact 
that the numerical results agree with the prediction by the 
three-layer model for the grain size of 0.01-1 |i.m, whereas 
the results agree with that by the two-layer model for the 
size > 10 [im. 

When there is scattering, some differences appear. For 
a — 0.01 [im (panel [a]), the scattering albedo u is negligi- 
ble in the wavelength interest (e.g., an effective temperature 
less than T* = 3,000 K in Figure 1). Thus, scattering vir- 
tually has no effect. For o = 0.1 u.m (panel [b]), ui for the 
stellar radiation is significant, but that for the diffuse radi- 
ation in the annulus (its effective temperature is less than 
about 300 K) is still negligible (see Figure 1). In this case, 
the temperature at the equatorial plane becomes slightly 
lowe r than that in the no s catteri ng case, which is consistent 
with lDullemond fc Nattal (|2003bh . For a = 1-10 urn (panels 
[c,d]), uj becomes significant for the radiation of the super- 
heated layer. In this case, we observe a plateau like structure 
at around Tp x * ~ 1 and a significant reduction of the equa- 
torial temperature. For a > 100 (j.m-1 mm (panels [e,f]), 
finally, u) becomes "grey" for all the radiation considered 
here. In this case, the temperature structure with scatter- 
ing becomes indistinguishable from that without scattering; 
a "grey" scattering has no effect on the temperature struc- 
ture in the optical depth coordinate. The physical reasons of 
these features will be discussed in section 3 with an analytic 
model. 

Even when the "grey" scattering, we find a difference 
of the temperature structures with/without scattering if we 
take the physical height as the coordinate as shown in Fig- 
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Figure 2. Vertical temperature structure in annuli with the radius of 1 AU but with various grain sizes, a, shown in each panel. The 
horizontal axis is the Planck mean extinction optical depth with the stellar effective temperature (T* = 3, 000 K). The solid lines are the 
no scattering case, the dotted lines are the case with the single scattering albedo for small wavelength u>o = 0.9, and the dashed lines are 
the case with u)q = 0.99. Two grey thick marks at the right-hand edge in each panel indicate the temperatures predicted by the analytic 
two-layer (upper mark; 93.3 K) or three-layer (lower mark; 78.5 K) models without scattering presented in section 3. 
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Figure 3. Same as Fig. 2, but the horizontal axis is the physical height from the equatorial plane normalised by the radius of the 
annulus. The diamonds indicate the lower boundary of the super-heated layer, at which the Planck mean extinction optical depth with 
the stellar effective temperature (T* = 3,000 K) is equal to the grazing angle (a = 0.05). 
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ure 3. The diamonds in Figure 3 indicate the lower boundary 
of the super-heated layer, at which the Planck mean ex- 
tinction optical depth with the stellar effective temperature 
(T* = 3, 000 K) is equal to the grazing angle. We find that 
with scattering, the height of the super-heated layer is al- 
ways enhanced; scattering causes more flaring disc (see also 
iDullemond fc Nattdl2003bh . This suggests that the scatter- 
ing may affect the global structure of the disc, which will be 
discussed in a future work. 



3 THREE-LAYER ANALYTIC MODEL WITH 
SCATTERING 

In order to understand the numerical results presented in the 
previous section, we here develop an analytic model as an 
extension of the seminal two-layer model by CG97: three- 
layer model with scattering. To describe fluxes across the 
boundaries of the layers, we adopt a two-stream Eddington 
approximation with isotropic scattering. The notations in 
this section are summarised in Table f . 



3.1 Model description 

Suppose two or three layers in an annulus as shown in Fig- 
ure 4. We assume that each layer is isothermal with the 
temperature determined by the radiation equilibrium in the 
layer. Then, we consider that each layer emits the radiation 
characterised by its temperature and other layers just work 
as absorption and scattering media for the radiation. The ra- 
diation from the central star is characterised by the stellar 
effective temperature. The characterisation of the radiation 
is done by the Planck mean and the characteristic frequen- 
cies are denoted by each subscript such as "*" for the stellar 
radiation (see the caption of Figure 4 and Table f). We de- 
note, for example, the extinction cross section characterised 
by the stellar effective temperature as k™ 4 . Note that all the 
quantities depending on the frequency are Planck averaged 
in this section. 

We call the two or three layers super-heated layer, mid- 
dle layer, and interior as shown in Figure 4. The super- 
heated layer is defined as only the layer exposed by the di- 
rect stellar radiation entering into the annulus with a small 
grazing angle a. The optical thickness of the layer is « a as 
shown by the numerical solutions in §2.2. Thus, we define 
the thickness of the super-heated layer as r s ext = K° xt E s = a, 
where E s is the gas column density of the layer. The inte- 
rior represents the isothermal part found in the numerical 
solutions. Thus, the boundary can be defined by the photo- 
sphere of its own radiation. However, we here simply define 
the interior as the part other than the super-heated and the 
middle layers. 

The middle layer is introduced by the following consid- 
eration. When the opacity coefficient decreases as the wave- 
length increases, the absorption of the radiation from the 
warm super-heated layer occurs well above the photosphere 
of the cold interior radiation. In this case, the interior is 
not warmed directly by the super-heated layer but by the 
"middle" layer where the radiation of the super-heated layer 
is effectively absorbed. We here define the thickness of the 



Table 2. Thickness of the upper two layers. 



middle layer as t, 1 



ext 
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Super-heated layer 
Middle layer 



_ext 
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^ext 
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Table 3. Condition of the two or three-layer models. 



Three-layer model 
Two-layer model 



-ext / t -CXt ^ 1 
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-ext /^.ext , , -i 
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arbitrary. On the other hand, when the opacity coefficient is 
grey, the middle layer with the above thickness is optically 
thick for its own radiation. Thus, the middle layer reaches 
the thermal equilibrium and is merged into the isothermal 
interior. In this case, we do not need to consider the middle 
layer. Therefore, we have two cases: the three-layer model 
< f and the two-layer model when 



f (or t; 



cxty 
fi/ m /-in 

ext 



f because the middle layer 



is merged into the interior). In other words, three layers 
are needed when f^'/C < 1 an d two layers are suffi- 
cient when ttf /k^ ~ f. Tables 2 and 3 are summaries 
of the thickness of the upper two layers and the condition 
of the two or three-layer models. Importantly, the seminal 
two-layer model is valid only when the opacity coefficient is 
grey in the frequencies interest. This fact has not seemed to 
be known well so far. 



3.1.1 Stellar fluxes 

When the grazing angle a is small, the stellar flux at the 
top of the annulus is 



Hi" = aW*B* 



(5) 



with the integrated Planck function B* = (o"sb/7t)T* and 
the dilution factor W, = fi*/47r, where Q* is the solid angle 
of the stellar photosphere from the top of the annulus. If 
only the fraction / v j s of the stellar photosphere is visible 
because of the optically thick disc, the solid angle becomes 
SI* = f v i B ir(R* I R) 2 , where R* is the stellar radius and R is 
the radius of the annulus. 

When there is scattering, a part of the incident stellar 
flux is reflected upwards a nd downwards by the super-heated 
layer. ICalvet et al.l (|l99fl ) presented an analytic expression 
of the scattered flux for i sotropic scattering. From equation 
(5) in lCalvet et all (|l99ll ), the outbound fluxes at the upper 
and lower boundaries of the super-heated layer (see Figure 4) 
become 



H? p = aW*B* 
and 



1 + X* 



1 + X. 



(6) 



(7) 



where u> 
quency and x* 



is the single scattering albedo at the stellar fre- 
In the derivation of equations (6) 
S s = a « 1, adopted 



and (7), we have assumed r s ex 



sity of the middle layer E m although this definition is rather a different upper boundary condition from ICalvet et al.l 
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Table 1. Notations in our analytic model. 



Notation Meaning 



Remarks 



a 

W* 
H7 

j^"down 
-^input 

T x 
B x 

^x 

t-abs 

K x 

^ext 

h x 

^-cxt 

U>x 

Xx 
r u P 
*x,y 
rdown 
1 x,y 

0>x,y 

bx,y 

Cx,y 

'I'input 

*i(2) 

*i(3) 



Grazing angle of the entering stellar radiation 
Solid angle of the stellar photosphere 

Dilution factor of the stellar radiation Q„/4tt 

Stellar input flux 

Upwards flux from the x layer 

Downwards flux from the x layer 

Total downwards flux from the super-heated layer 

Temperature of the x layer 

Frequency integrated Planck function with T x (ag-Q/ir)T^ 
Gas mass column density of the x layer 

Absorption cross section per unit gas mass for the radiation with T x 
Extinction cross section per unit gas mass for the radiation with T x 
Extinction optical depth of the x layer for the radiation with T y 
Scattering albedo for the radiation with T x 
Square-root of the thermal coefficient 1 — ui x 
Upwards intensity of the radiation with T y from the x layer 
Downwards intensity of the radiation with T y from the x layer 
Thermal coefficient in the x layer for the radiation with T y eq 
Reflection coefficient in the x layer for the radiation with T y eq 
Transmission coefficient in the x layer for the radiation with T y eq 
Ratio of Zfjnput to H'™ eq 
Reduction factor of B\ by scattering in the two-layer model eq. 
Reduction factor of B\ by scattering in the three-layer model eq 



^cxt V 

K y ^ x 



(B13) 

(B14) 

(B15) 

(20) 

(23) 

(29) 



Subscripts x and y are * for stellar quantities, s for super-heated layer quantities, m for middle layer quantities, or i for interior 

quantities. 
crsB i s the Stcfan-Boltzmann constant. 
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Figure 4. Schematic diagram of the two-layer and the three-layer models. We are considering an annulus clipped from a flaring disc 
surrounding a young star as shown in the left-hand picture, fn the annulus, we consider two or three layers as shown in the right-hand 
pictures. Each layer is assumed to be isothermal with a temperature T. The fluxes crossing the boundaries of layers are denoted as H up 
or H dovJn depending on their direction . We set the direction of arrows in the right-hand pictures as the positive direction for the fluxes. 
The subscript of each quantity indicates the layer associated with the quantity: "s" for the super-heated layer, "m" for the middle layer, 
and "i" for the interior. The quantities with the subscript of "*" arc associated with the radiation from the central star and H'™ is the 
stellar flux at the top of the annulus. See also Table 1 for notations. 
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and neglected the term e 1 for the downwards 
flux. Note that the total scattered flux is 77, p + 77 d own = 
w,aff,B, = Lj*Hi n and i?" p = 77 down = if w* = (no 
scattering case). 

3.1.2 Super-heated layer fluxes 

The radiation characterised by the temperature of the super- 
heated layer T a is produced only in the super-heated layer. 
In other layers, this radiation is not produced but is ab- 
sorbed or scattered. The super-heated layer is vertically 
optically thin for its own radiation as r s cxt = Kg Xt E a = 
( K s Xt / K * xt )o ! •C 1 because the grazing angle a is small and 
we have (Kf7 K ?') ^ 1- On the other hand, the total opti- 
cal depth of other layers is very large. Thus, we consider a 
geometry that a thin isothermal layer lies on a semi-infinite 
absorption and scattering slab. 

The above of the super-heated layer is assumed to be 
vacuum; there is no downwards input radiation with the 
temperature T s at the top of the layer. However, there is 
upwards input radiation at the bottom of the layer because 
of the reflection by the semi-infinite interior below the layer. 
The upwards and downwards radiation intensities from the 
super-heated layer (7^ P and 7 d ° wn , respectively) in the two- 
stream Eddington approximation (the cosine of the angle 
between the stream lines and the normal of the layer is set 
to be ±1/V3) become from equations (B12) and (B16) 

7 s up = a s , s B s + Cs tS I^ , (8) 



and 



= a s ,sB E + b B , B I" 



(9) 



where 7 ; up is the upwards input intensity from the interior, 



and a s 



and c s s are the thermal, reflection, and trans- 



mission coefficients in the super-heated layer for the radi 
ation characterised by the temperature T s (see eqs. [B13- 
B15]). The reflected intensity from the interior becomes 

rdown 



(10) 



where bi, s is the reflection coefficient in the interior for the 
radiation with T s . Note that the interior itself does not emit 
radiation with T s . 

We can always obtain the unique exact solution of 
^ P , 4 d ° wn , and 7" s p from equations (8-10), whereas we de- 
rive an approximate solution of them in this paper. As 
found equation (B18), for a semi-infinte slab, we have 
6i,s ~ (1 — Xs)/(1 + Xs), where Xs = VI — u s with uj s be- 
ing the scattering albedo for the radiation with T s . Since 
the super-heated layer is optically thin for its own radia- 
tion, we can approximate 6 S , S <C 1 and c SjS « 1. We also 
find a s , s » V3Xs T sT = V^Xs ( K s Xt / K * xt ) a from equation 
(B13) for a small optical depth. Then, we obtain 7" P ss 
(1 + b Us )a s , s B s , 7 down « a s , s B s , and I** « bi, s a s , s B s . The 
upwards and downwards fluxes from the super-heated layer 
are 77 s up = 7^/(2^3) and 77 down = (7*°™ - 7^/(2^3). 
Therefore, we obtain 



77 s up = aB s 
and 



77 s = aB s 





[ Xs 2 1 




.1+Xs. 



X., 



1+Xs 



(11) 



(12) 



3.1.3 Middle layer fluxes 

As discussed in the section 3.1, we consider the middle layer 
only when the layer is optically thin for its own radiation. 
Although the middle layer is sandwiched between the super- 
heated layer and the interior, we neglect the effect of the 
super-heated layer because the optical depth of the layer is 
very small (see above) . The optical depth of the interior is so 
large that we can regard it as a semi-infinite medium. Thus, 
we have the same setting as the super-heated layer, other 
than the optical thickness of the layer, r™m = (C t /«f It ). 
Following the section 3.1.2, we have 



Hit? = Br, 



and 



rrdoWll t-> 



Xn 



1+Xi 





Xm 


V Ki xt y 


1 + Xm_ 



(13) 



(14) 



where 73 m is the integrated Planck function with the tem- 



perature of the middle layer T m 



is the scattering albedo 

VI - Wra- 



for the radiation of the middle layer, and \ 



3.1.4 Interior flux 

We always consider the interior to be optically thick for 
its own radiation. On the other hand, other layers above 
the interior are considered to be always optically thin. If 
we neglect the effect of the upper layers, the interior is re- 
garded as a semi-infinite isothermal medium without inci- 
dent flux of the radiation with the temperature T|. In this 
case, based on equation (B12), the upwards outbound flux 
becomes 77 ; up = ay.Bi/(2-\/3) with a\ i being the thermal co- 
efficient in the interior for the radiation with T; and 7?i being 
the integrated Planck function with 71. For a semi-infinite 
medium, ay ~ 2xi/(l + xO; where Xi = VI — w i with u>i 
being the single scattering albedo for the interior radiation 
(eq. [B17]). Therefore, we have 



77" 



V3 



Xi 



1 + Xi 



(15) 



Note that the factor in [ ] is equal or less than 1/2. The 
equal is true without scattering. Therefore, in general, scat- 
tering reduces the radiation energy loss from the interior, 
i.e. the green-house effect. 



1 ICalvet et all 1,199 1[) adopted J(0) = 277(0), where J(0) and 
7f(0) are the mean intensity and the mean flux of the scattered 
stellar radiation at the top of the medium. On the other hand, 
we adopted J(0) = V377(0). This difference is due to the choice 
of the angle of the stream line. 



3.2 Temperature of the super- heated layer 

The radiation energy conservation in the super-heated layer 
is 



77i n - 77^' 



rrdown 



77 s u 



77 s a 



(16) 
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Figure 5. Temperature of the super-heated layer as a function 
of grain size: (a) annulus with the radius of 1 AU, (b) 10 AU, 
and (c) 100 AU. The numerical solutions are shown by symbols; 
circles are no scattering case, triangles are the case with the single 
scattering albedo for small wavelength u)q = 0.9, and squares are 
the case with u>o = 0.99. The solid lines are the analytic model 
given by equation (18). 



The left-hand side is the net input flux of the stellar radia- 
tion and becomes aW,B,(l-u,). When there is scattering, 
the input flux is reduced by a factor of 1 — w*. The right- 
hand side is the net output flux of the super-heated layer 
and becomes a_B s (Ka Xt /K* xt )(l — u s ). Thus we obtain 



B s = B,W* 



K% hs 



1/4 



1/4 



(17) 



(18) 



where we have applied K abs = (1 — cj)k cx *. Therefore, we ex- 
pect T s to be independent of the existence of scattering. This 
result can be recovered by a microscopic consideration for 
each dust grain in the radiation equilibrium with the stellar 
radiation field in the super-heated layer. We also expect that 
T s shows two asymptotic values since (Kj bs /reJ bs ) = T*/T s 
for a small grain size and («;* bs //ts bs ) = constant for a large 
grain size as shown in Figure 1. 

Figure 5 shows the temperature of the super-heated 
layer as a function of grain size. The numerical solutions pre- 
sented in section 2 are shown by symbols: circles for the no 



scattering case, triangles for the u>o = 0.9 case, and squares 
for the too = 0.99 case. The analytic model presented in 
equation (18) is shown by the solid lines. To obtain the 
solution T s of equation (18), we need an iterative proce- 
dure because Kg hB depends on T B . Figure 5 first shows that 
the scattering hardly affect T s as expected by the analytic 
model. Indeed, the cases of the three different albedos are su- 
perposed for almost all grain sizes, although a temperature 
enhancement of a few percent is observed in the scattering 
cases at a grain size of 0.1-1 um. Figure 5 second shows 
that the numerical solutions can be divided into two cases: 
higher temperature for smaller grain size (< 0.1 p) and 
lower temperature for larger grain size (> 10 um), which is 
also expected by the analytic model. The agreement between 
the numerical solutions and the analytic model is excellent 
although the numerical solutions give about 3% higher tem- 
perature than the analytic model. This small difference is 
probably caused by neglecting the absorption of the interior 
radiation in the super-heated layer in the analytic model. 



3.3 Temperature of the interior 

We can define the radiation flux input into the interior (in- 
cluding the middle layer) of the annulus as 



//input 

where 



H t +H S = aW4B*$m P ut 



nput 



+ 



(19) 



(20) 



and we have eliminated B s by equation (17). If there is no 
scattering (w* = lj b = 0, i.e. \* = Xs — 1), we have $i nput = 
1/2. That is, f2"i Ilput = aW,B,/2 = H'T/2; a half of the radi- 
ation energy received at the top of the a nnulus is input into 
the interior l|Chiang fc Goldreichl [l997h . For the isotropic 
scattering case, we always have //input < Hl n /2; isotropic 
scattering always reduces the energy input into the inte- 
rior re lative to the no scattering case (|Dullemond fc Nattal 
l2003bh . 



3.3.1 Two-layer model 

The energy conservation in the interior under the two-layer 
model is 



//input — H\ p . 

From equations (15) and (19), we obtain 

Bi = VSaW,B^ i(2 ) , 

where 

'l + Xi 



4>i 



(2) 



Xi 



nput • 



The interior temperature becomes 
T = T, [\/3aW/,$i (2) ] 1/4 . 



(21) 
(22) 

(23) 
(24) 



Therefore, T is proportional to a factor of ■ When there 
is no scattering (c&input = 1/2 and \i — 1), we have <& i ( 2 ) = 1. 

As discussed in section 3.1, the two-layer model is valid 
if the opacity coefficient is regarded as grey at all the fre- 
quencies interest: Kf xt /«s Xt ~ 1. If this condition is satisfied, 



Vertical temperature structure of protoplanetary discs 9 



the scattering is also grey; o>» = cu s = ui\ (i.e. x* — Xs = Xi)- 
In this case, we have 

$i( 2 ) = 1 (no/grey scattering) . (25) 

Importantly, the factor $i(2) for grey isotropic scattering is 
independent of the albedo and equal to the case without 
scattering. This is caused by the fact that the reduction 
of the flux input into the interior by scattering ($i npu t in 
eq.[19]) is completely offset by the reduction of the flux out- 
bound from the interior by scattering ((1 + Xi)/Xi m ec l-[15]) 
for grey and isotropic scattering. Therefore, we expect the 
same interior temperature for grey and isotropic scattering 
case as that for no scattering case in the two-layer model. 



3.3.2 Three-layer model 

The energy conservations in the middle layer and the interior 
under the three-layer model are 

-ffinput = #m P + i/^° Wn , (26) 

and 

= #j up . (27) 

Since #£ own = XmC from equations (13) and (14), we 
obtain H? p = X m/(1 + Xn,)H input . Thus, 

Bi = v / 3aW»£,$i( 3 ) , (28) 
where 

'?Cm\ ( 1 + Xi 



<E>i(3) = 



4>i 



nput ■ 



Xi / V 1 + Xn 

The interior temperature becomes 
Ti = T* [V3aW*$ m ] 1/4 . 



(29) 



(30) 



.1/4 



Therefore, Ti is proportional to a factor of $J 3 y 

From the discussion in section 3.1, the three layer model 
is valid when re™* / re^ < 1. This condition is satisfied when 
the grain size, a, is smaller than 10-100 |xm. When a < 0.01 
(i.m, the scattering albedo is negligible at the all frequencies 
interest. In this case or the no scattering case {uj x = 0, i.e. 
Xx = 1)> we have 
1 
2 

Comparing this with the two-layer model, we find that Ti 
in the three-layer model is reduced by a factor of (1/2) 1 / 4 
relative to that in the two-layer model even for no scattering 
case. 

When 0.01 < a < 0.1 urn, loi w oj m w lu b « (xi ~ 
Xm ~ Xs ~ 1) but w* > (x* < 1). In this case, we have 

X* 



<E>i(3) = - (no scattering) . 



(31) 



*i(3) = fC2~X*)- 



(32) 



Thus, we find $1(3) < 1/2; we expect a reduction of Ti rela- 
tive to that without scattering. When 0.1 < a < 1-10 [im, 
Wi » w m « (xi ~ Xm ~ 1) and w s w w* > (xs ~ X* < !)■ 
In this case, we have 



(3) 



X* 



i + x. 



(33) 



When 1-10 < a < 10-100 ^.m, lj { w (xi ~ 1) and w n 
w s » w» > (xm ~ Xs ~ X* < !)• Then, we have 
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Figure 6. Reduction factor of the interior flux by isotropic scat- 
tering in the analytic three-layer model, $1,(3), as a function of 
the albedo for the stellar radiation, lu*. The solid, dot-dashed, 
and dashed lines indicate the cases with a typical grain size of 
°- 01 £ a £ 0- 1 0- 1 £ a £ i- 10 M™, and 1-10 < a < 10-100 
|jm, respectively. The reduction of the interior temperature rel- 
ative to that in the two-layer model without scattering is given 
1/4 

by <f>. '.g. . Note that the interior temperature in the three-layer 

model is reduced by a factor of (1/2) 1 / 4 even for (j* = (i.e. no 
scattering) . 



4>i 



(3) 



1+X* 



(34) 



When a > 10-100 (i.m, the opacity becomes almost grey, 
thus, the three-layer model is no longer valid. We should 
choose the two-layer model in this case. 

Figure 6 shows the scattering reduction factor $i,(3) 
given by equations (32-34) as a function of the albedo for 
the stellar radiation. We find that the factor decreases from 
equation (32) to (34), in other words, as a function of the 
grain size. The factor gives the reduction of the inte- 

rior temperature, Ti, in the three-layer model by isotropic 
scattering relative to that in the two-layer model without 
scattering. For typical grain sizes found in protoplanetary 
discs of 0.1-10 (im, equation (33) would give a good approx- 
imation for the reduction factor. If a;* w 1, equations (32) 
and (33) are reduced to w x* = \/l — (J*. In this case, we 
expect the reduction factor of Tj to be ~ (1 — o;,) 1 ^ 8 which 
is found in Figure 7 later. 



3.3.3 Comparison between numerical and analytic results 

Figure 7 shows the temperature of the interior as a func- 
tion of the grain size. The numerical solutions presented in 
section 2 are shown by symbols: circles for the no scatter- 
ing case, triangles for the coo = 0.9 case, and squares for 
the wo = 0.99 case. The analytic models are shown by the 
solid lines. In analytic models, we always assume the opti- 
cally thick interior. Most of the numerical solutions shown in 
Figure 7 are really optically thick. However, when the grain 
size is larger than about 100 [im, because of the reduction 
of «o bs given by equation (4), some cases without scattering 
and with ojq — 0.9 become optically thin. In such cases, we 
find relatively high temperatures. 

The solutions of the analytic models are obtained as 
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Figure 7. Temperature of the disc interior expected in isotropic 
scattering case as a function of grain size: (a) annulus with the 
radius of 1 AU, (b) 10 AU, and (c) 100 AU. The numerical so- 
lutions are shown by symbols; circles are no scattering case, tri- 
angles are the case with the single scattering albedo for small 
wavelength ujq = 0.9, and squares are the case with u>o = 0.99. 
The solid lines are the analytic model described in sections 3.3.1 
(two-layer model; right-hand side) and 3.3.2 (three-layer model; 
left-hand side). We connect these two models around a grain size 
indicated by the dotted line, where a jump appears because of this 
connection (see text in detail). Upwards deviations of numerical 
solutions from the two-layer model found in panels (b) and (c) 
are caused by that the interiors of these cases are optically thin 
for these own radiation. 



follows: for the two-layer model, we solved equation (24) to 
obtain 71. We need an iterative procedure because the term 
<3>i(2) depends on 71. For the three- layer model, we obtained 
71 from equation (30) after obtaining T m from 



aW*B* 



$input 



1/4 



(35) 



which is derived from equations (13), (14), and (26). We 
again need an iterative procedure to obtain both of T m and 
71. 

As discussed in §3.1, the two-layer model is valid only 
when «:f xt //tg Xt « 1. Otherwise, we should adopt the three- 
layer model. Here, we connect these two models at the grain 
size where K° xt /Kg Xt = 0.8 in the two-layer model, for ex- 
ample; we adopt the two-layer model for a larger grain size 
than it and adopt the three-layer model for a smaller grain 



size than it. This threshold is rather arbitrary, but we find 
that this choice is good as seen in Figure 7 where the con- 
nected analytic models reproduce the numerical solutions 
reasonably well. Note that a sudden jump on the solid lines 
in Figure 7 is caused by this connection and numerical so- 
lutions also show relatively rapid change of the temperature 
around there. 

When there is no scattering, we find two asymptotic 
values of 71. Interestingly, T\ for smaller grain size is lower 
than that for larger grain size, whereas T s show the opposite 
trend in Figure 5. Note that the radiation flux input from 
the super-heated layer is in fact independent of T s as shown 
in equation (20) in the no scattering case. Nevertheless, the 
numerical solutions show a factor of (1/2) 1 / 4 reduction of 
71 for small grain cases. This is excellently explained by the 
three-layer model as found in equation (31). The physical 
explanation is as follows: when the grain size is small, the 
radiation flux from the super-heated layer is absorbed by 
the middle layer once. Then, the middle layer emits a half 
of the absorbed energy downwards but the rest of the half 
goes upwards. Therefore, the interior does not receive a half 
of the stellar radiation energy absorbed by the super-heated 
layer but receive only a quarter of the energy when the grain 
size is small. 

When there is isotropic scattering, we find further re- 
duction of 71 in a range of the grain size a = 0.1-10 or 100 
[im. For a ~ 0.01 u.m, the scattering effect is not observed 
because of negligible albedo. For a > 10-100 [im, the scat- 
tering effect is not observed, either (but except for the op- 
tically thin cases). This is nicely explained by the two-layer 
model with grey isotropic scattering as found in equation 
(25). The physical reason is that the grey isotropic scat- 
tering equally reduces the downwards flux from the super- 
heated layer (eq. [19]) and the upwards flux from the interior 
(eq. [15]). The reduction of 71 due to scattering for a — 0.1- 
10 vim is explained by the three-layer model as summarised 
in equations (32), (33), and (34). The physical reason of the 
71 reduction is that the scattering reduces only the down- 
wards fluxes from the super-heated layer and from the mid- 
dle layer because the scattering albedo at the frequency of 
the interior radiation is still negligible. 



4 CONCLUSION 

We have examined the effect of scattering of the diffuse ra- 
diation on the vertical temperature structure of protoplan- 
etary discs. This is motivated by the fact that scattering 
albedo increases as the size of dust grains grows in the discs. 
In particular, large icy grains have a significant albedo even 
in the infrared wavelength. For this aim, we have developed a 
ID plane-parallel numerical radiation transfer code includ- 
ing isotropic scattering of the diffuse radiation as well as 
that of the incident radiation. We have also developed an 
analytic model with isotropic scattering both of the diffuse 
and the incident radiations in order to interpret the solu- 
tions obtained from the numerical simulation. All results of 
the numerical simulations has been nicely reproduced by the 
analytic model. 

The analytic model presented in t his paper is an exten- 
sion o f the seminal two-layer model bv lChiang fc Goldreichl 
l| 19971 ); we have introduced a new layer between the super- 
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heated surface layer and the disc interior. This middle layer 
(or disc atmosphere) is required when the absorption of the 
radiation from the super-heated layer occurs well above the 
photosphere of the opaque isothermal interior. This situa- 
tion is realised if the dust opacity is negatively proportional 
to the wavelength. Thus, we should consider three layers 
rather than two layers if the grain size is smaller than about 
10 u.m. On the other hand, for grey opacity, which is realised 
if the grain size is > 10 urn, the standard treatment with 
two layers is justified. 

We have found from the numerical simulation that the 
dust temperature of the disc surface is almost not affected by 
scattering. This is because the temperature is determined by 
the radiation equilibrium of grains in the incident radiation 
field locally. The grain size has an effect on the surface tem- 
perature via the wavelength dependence of the dust opacity. 
There are two asymptotic temperatures: the higher temper- 
ature is realised by small grains (< 0.1 (i.m) and the lower 
temperature is realised by large grains (> 1-10 [im). This 
is because the emission efficiency relative to absorption ef- 
ficiency of the small grains is smaller than that of the large 
grains which have grey opacity. This trend of the numerical 
solutions is excellently reproduced by the analytic model of 
the super-heated layer. 

The numerical simulations without scattering show that 
the dust temperature of the optically thick interior also 
has two asymptotic values: the lower temperature for small 
grains (< 0.1 |i.m) and the higher temperature for large 
grains (> 10-100 (i.m). Thus, the trend is opposite from 
the surface temperature. The higher asymptotic tempera- 
ture for large grains is well matched with the prediction of 
the two-layer model. On the other hand, the lower asymp- 
totic temperature for small grains is a factor of (1/2) 1//4 
lower than the prediction of the two-layer model. In fact, 
the flux input from the super-heated layer is always same al- 
though the interior temperature is different as a function of 
the grain size. This phenom enon has been already found by 
iDullemond fc Nattal l|2003al ) who attributed it to the energy 
loss fro m the interior by the rad iation at a long wavelength 
(see also IDullemond et al.ll2002r ). We have proposed another 
interpretation by the middle layer between the super-heated 
layer and the interior; the super-heated layer gives a half of 
the absorbed energy of the incident radiation to the mid- 
dle layer which gives a half of the obtained energy to the 
interior. This three-layer model exactly predicts a factor of 
(1/2) 1 / 4 reduction of the interior temperature for the small 
grain case. 

The numerical simulations with isotropic scattering also 
show two asymptotic temperatures of the interior. Interest- 
ingly, these asymptotic temperatures of no scattering and 
isotropic scattering cases are the same. For small grains 
(< 0.1 nm), since the scattering albedo is negligible for all 
wavelengths interest, the same temperature is trivial. The 
same temperature for large grains (> 10-100 \im) has been 
nicely explained by the two-layer model with grey opacity. 
The physical mechanism is the exact offset between the re- 
duction of the flux input into the interior by scattering in the 
super-heated layer and the reduction of the flux output from 
the interior by scattering in itself (i.e. green-house effect). 

For grain sizes of 0.1-10 urn, which are expected by a 
moderate growth of the size in the discs, we have found a fur- 
ther reduction of the interior temperature in isotropic scat- 



tering cases relative to that without scattering from the nu- 
merical simulations. This reduction has been well explained 
by the three-layer analytic model. The physical mechanism 
is the wavelength dependence of albedo; the flux input into 
the interior is reduced by scattering in the super-heated and 
middle layers, whereas the flux output from the interior is 
not reduced because of negligible (or weak) scattering. Note 
that the interior flux has a longer wavelength typically, thus, 
the albedo for the radiation is smaller. 

In conclusion, the scattering of the diffuse radiation can 
affect the vertical temperature structure of protoplanetary 
discs significantly when the grain size grows to be about 1— 
10 (xm. We need to investigate the effect in a global disc 
model in future. The analytic model presented in this paper 
could be useful to understand the physics determining the 
temperature structure in the discs. 
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APPENDIX A: VARIABLE EDDINGTON 
FACTOR METHOD WITH ISOTROPIC 
SCATTERING IN PLANE-PARALLEL SLAB 

We here present our numerical radiation transfer method in 
a plane-parallel m edium in detail. The me thod is based on 
that developed bv lDullemond et al.l (|2002T >. but is extended 
to treat isotropic scatterings of both of the incident radiation 
(from the central star) and the diffuse radiation (from dust 
grains in the medium). We find a solution, in which the 
radiation field, the temperature structure, and the density 
structure are consistent with each other, iteratively as the 
following procedure: 

(i) assuming an initial temperature structure 

(ii) solving the density structure consistent with the given 
temperature structure under the hydrostatic equilibrium 

(iii) solving the transfer of the incident (stellar) radiation 
with the grazing angle recipe 

(iv) solving the transfer of the diffuse radiation with a 
variable Eddington factor method and obtaining the tem- 
perature structure under the radiation equilibrium 

(v) checking the convergence of the temperature structure 
and if not going back to the step (ii) 

In the following we describe the set of equations and as- 
sumptions and the result of a benchmark test. 

Al Hydrostatic equilibrium 

Suppose an annulus clipped from a protoplanetary disk to be 
a plane-parallel medium. We set the coordinate z as the ver- 
tical height of the medium. The origin z — is the equatorial 
plane of the annulus and we set a mirror boundary condi- 
tion there. Assuming the vertical hydrostatic equilibrium, 
we obtain the density of gas in the medium p(z) consistent 
with the temperature structure T(z) which is assumed as 
an initial guess or is obtained by the previous step of the 
iteration. We assume that the gas temperature is the same 
as the dust temperature which is determined by the radia- 
tion equilibrium. This assumption is usually well established 
in the protoplanetary disc because the collision between gas 
particles and dust grains occurs enough frequently. 
The vertical hydrostatic equilibrium is given by 

= (AD 

where P is the gas pressure and g is the gravitational ac- 
celeration. In a protoplanetary disc, the self-gravity of the 
disc is negligible relative to the gravity of the central star. 
Thus, we have g = GM^z/ R 3 , where G is the gravitational 
constant, M* is the mass of the central star, and R is the 
distance from the star (or radius of the annulus considered). 



We have assumed R^> z. The gas pressure is given by the 
equation of state for the ideal gas as P — (pfcBT)/(/i m m p ), 
where /x m is the mean molecular weight, m p is the proton 
mass, and &b is the Boltzmann constant. Then, equation 
(Al) is reduced to 



din p 



z + 



dlnT 



/ fj,m p \ / GM t \ 
dz ' \k B TJ \ R 3 J ~ </.: 



(A2) 



If we integrate equation (A2) from z — with a given T(z), 
we obtain the functional shape of p{z). The absolute value 
of p(z) is scaled by 



E = 2 



p(z)dz, 



(A3) 



where E is the gas column density and z max is the maximum 
height for the numerical calculation. We set z max = R in our 
calculation and we set the minimum value of p — 10 -25 g 
cm -3 just for avoiding a too small value of the density. 

We set the optical depth coordinate r for the radiation 
transfer from the obtained p(z) as 



Tv («) = 



p{z)K^ t dz , 



(A4) 



where rej xt is the extinction cross section by dust grains per 
unit gas mass at the frequency v. 

In fact, we can obtain the temperature structure as a 
function of the optical depth T(r) by the radiation transfer 
without the density structure as a function of the vertical 
height p(z). The reason why we calculate p(z) and the rela- 
tion between the optical depth and the vertical height t(z) 
is to see T(z) as shown in Figure 3. Another reason, which 
may be more important, is to determine the grazing angle 
consistent with the disc global structure for future calcula- 
tions. 

A2 Transfer of the incident radiation 

Let us consider an incident radiation beam entering the 
plane-parallel medium. The angle between the incident ray 
and the surface of the medium is called the grazing angle, 
a. In a protoplanetary disc, the grazing angle a is usually as 
small as ~ 0.05 radian fe.g.. |P'Alessio et al.| [2006 V Thus, we 
use the approximation sin a ~ a. The optical depth along 
the incident ray becomes T v (z)/a. Thus, the mean intensity 
of the (direct) incident radiation is given by 



Jt (z) = J* 



(A5) 



where J* max is the mean intensity at the top of the medium 
(i.e. z = z max ). For the incident radiation from the central 
star, j; max = B„(T,)n,/4n, where B„{T,) is the Planck 
function with the stellar effective temperature T* and f2* = 
7r(i?*/7?) 2 / v is is the solid angle of the stellar photosphere 
visible from the top of the medium (/ v i s is the visible fraction 
of the stellar photosphere). Finally, we give the absorbed and 
extincted energy density of the incident radiation per unit 
time interval at the height z as 



gabs (2) = 



and 



p(z)/v abs 47r dv , 



g ox t(z) = / p{z)K^ t A-Kj'*dv , 



(A6) 



(A7) 
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where K^ bs an d «^ xt are the absorption and the extinction 
cross section by dust per unit gas mass at the frequency v. 

A3 Transfer of the diffuse radiation 

The transfer equation of the diffuse radiation (or the radia- 
tion reprocessed by dust) in a plane-parallel medium is 

dluu. 



dz 



cxt t i ext £i 

-pK v li/fj, -r pt^ir ; 



(A8) 



where p is the cosine of the angle between the ray and the 
z coordinate, I v/1 is the specific intensity at the frequency v 
towards the direction p, and S„ is the source function which 
is given by 



Sv = (1 — UJu)B v (T) + UJ V J V + UJvJl 



(A9) 



where u v is the single scattering albedo at the frequency 
i>, B V (T) is the Planck function with the dust temper- 
ature T, and J„ is the mean intensity of that is, 

1 f 1 

J v = — I Iv^dfi. We have assumed that the scattering is 

2 J-i 

isotropic for the simplicity. Note that we consider the scat- 
tering of the diffuse radiation as the second term in equation 
(A9). The third term accounts for the scattering of the in- 
cident radiation. To determine the dust temperature T, we 
assume the radiation equilibrium as 



pC s B v {T)du = 



abs t j , Qabs 

pK v J v dv + — . 

47T 



(A10) 



To obtain the mean intensity J„, we adopt a variable 
Eddington factor method. The first and second moments of 
equation (A8) are 



dH v ab 

—j— = pn v (B v — J v ) + pn v J u . 

and 

dK v 



dz 

where H v 



cxt j j 

= -pn v H v , 



(All) 



(A12) 



Iv^pdji and K v 



Iv^p dp,, and we 



have used equation (A9) to eliminate S v . Note that K^ hB = 



(1 



and the scattering cross section per unit gas 



mass Kj ca = ujvK^ 



If we integrate equations (All) and 
(A12) over the frequency, we obtain 



dH 

dz 
and 

dK 

dz 



Qcxt 

4tt 



pn^Hjjdu , 



(A13) 
(A14) 

K v du, and we have 



where H = H v dv and K = 

Jo Jo 
used the radiation equilibrium (eq.[A10]) and equations (A6) 

and (A7) in equation (A13). Finally, we introduce the Ed- 
dington factor as the closure equation: 



h = 



K 



J ' 
where J = 



(A15) 



J v dv. 



200 
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100 - 
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Figure Al. Result of a benchmark test without scattering pro- 
posed by C. P. Dullemond. The circles are the reference solutions 
by C. P. Dullemond, whereas the solid line is our solution. 



We do not assume /e = 1/3 (constant) as in the 
usual Eddington approximation, but obtain /e directly from 
I vll which is calculated by the formal solution of equation 
(A8). The integration of the formal solution is performed 
with a parabolic i nterpolation of S v am ong three succes- 
sive spatial point s llOlson &: Kunasj 19871 ) . As described in 
iDullemond et all l|2002f ). we alternate the integration of the 
formal solution (i.e. ray-tracing) with the integration of the 
moment equations (eqs. [A13]-[A15]) and the determination 
of T by equation [A10] until we reach a convergence in T (the 
difference between the two successive iterations becomes less 
th an . 1%). In addition, we adopt the acceleration algorithm 
byE3|l9Zi) for a rapid convergence. For no scattering case, 
the convergence is very rapid independent of the total op- 
tical depth of the medium, typically 20 iterations. On the 
other hand, a factor of ~ 10 times more iterations depending 
on the albedo are needed for isotropic scattering case. 



A4 Benchmark test 

Figure Al shows the result of the benchmark test No. 3 pro- 
posed by C. P. Dullemond on his web pageQ. The settings 
are as follows: the stellar temperature T„ = 3, 000 K, the 
stellar radius R, — 2.0 Rq, the distance from the star R — 1 
AU, the grazing angle a = 0.05, the visible fraction of the 
stellar photosphere / v ; s = 0.5, the total visual optical depth 
of the disc To.55(j.m = 10 4 (i.e. ro.ssum = 5x 10 3 at the equa- 
torial plane), and no scattering. We used the dust opacity 
model downloaded from the web page which is the same as 
that assumed in the calculation by C. P. Dullemond. The 
figure shows an excellent agreement between both results. 



http: / /www. mpia-hd.mpg.de/homes/dullemon/radtrans/benchmarks/ 
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Figure Bl. Plane-parallel isothermal medium with isotropic 
scattering. 

APPENDIX B: ANALYTIC EXPRESSION OF 
RADIATION FIELD IN ISOTHERMAL, 
ABSORPTION, AND ISOTROPIC 
SCATTERING MEDIUM 

Here, we derive an analytic expression of the radiation field 
in an isothermal medium with absorption and isotropic scat- 
tering. In the derivation, we adopt th e Eddington approx- 
imati on with two-stream lines (e.g., iRvbicki fc Lightmanl 
1 19791 ). We consider a diffuse incident radiation. 

Suppose a plane-parallel medium (see Figure Bl). We 
set the extinction (absorption+scattering) optical depth co- 
ordinate r along the normal of the medium. The total ex- 
tinction optical depth of the medium is set to be tl. Then, 
suppose that the single scattering albedo in this medium is 
to and this medium is isothermal and in the thermal equilib- 
rium. The thermal radiation is denoted as B. Let us consider 
two-stream lines with the direction of /i — ±l/\/3, where /x 
is the cosine of the angle between the ray and the r coor- 
dinate. When the direction of the rays is denoted by the 
subscript of + or — , the equation of the radiation transfer 
becomes 

1 dl± 



— - , = S± — I± . 

The source function can be expressed as 
S± = {l-u)B + uJ, 



(Bl) 



(B2) 



where J is the mean intensity. In the two-stream approxi- 
mation, we can define the mean intensity J, the mean flux 
H, and the mean radiation pressure K as 



J=\{I+ + I-), 
and 

K= 1 -(I + + I.)= 1 -J. 

The first and second moments of (Bl) with the source 
function (B2) are 

dH 

~dT 

and 



J), 



(B3) 
(B4) 

(B5) 



(B6) 
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Figure B2. Thermal, reflection, and transmission coefficients as 
a function of effective optical depth of an isothermal, absorption, 
and isotropic scattering medium. In each panel, we show three 
cases of the single scattering albedo: u) = (solid curve), 0.90 
(dotted curve), and 0.99 (dashed curve). The reflection coefficient 
b is always zero when ut = 0, so that we cannot see the solid curve 
in the panel (b). The effective optical depth = y/ 3(1 — tl 
if the total (absorption+scattering) optical depth of the medium 
is tl . Asymptotic values of the coefficients are given in equations 
(B17)-(B19). 



ldJ 

3 dr 



= -H. 



(B7) 



With an effective optical depth r c = r-y/3(l — lj), we have 
d 2 J 



= J-B, 



(B8) 



from equations (B6) and (B7). 

If we have incident radiations at the upper and lower 
boundaries as Iq 1 and 7 L n , respectively, the boundary condi- 
tions are 



j ( o) = # + vr^(£) o , 

and 

j(T C L ) = /["-yr-^(£) L , 



(B9) 



(B10) 



where = tlv^3(1 — u>). The solution of equation (B8) 
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with the boundary conditions (B9) and (BIO) is 



J(t.) = S 



1 - 



1 + X + (1-X)e- T ^ 
Jin (l + X)e- T --(1-XK 2T ° L+T ° 
° (1 + X ) 2 -(1-X) 2 e- 2 ^ 

+ g (1 + x)e " T ° +T °- (1 ' x)e f^, (Bll) 

(1+ X )2-(1- X )2 e -^ 



where x = vl — w. In the case without incident radiation 
(i.e. /g n = J?" = 0)i this solution is ex actly same as equa- 
tion (28) of iMivake fc Nakagawal i|l993| y 

Then, let us consider the outbound intensity at the sur- 
face of the medium. At r = 0, the outbound intensity can 
be 

I out = aB + bll n + cl'l n , (B12) 

where a, b, and c can be called as "thermal", "reflection", 
and "transmission" coefficients, respectively. Since J(0) = 
(I/2)(4 n +/o ut ), we obtain from equation (Bll) 



2x(l 



1 + X + (1 - X)e- T " 



(B13) 



(i-x)(i + x)(i- 



e 



-2rt 



b = ^ A/v A/v =4 , (B14) 

(1 + X) 2 -(1-X) 2 e- 2 ^ L 

and 



(B15) 



(1 + X) 2 - (1 - X) 2 e- 2 ^ 
At r = tl, we obtain symmetrically 

J~* = aB + bll n + c4 n . (B16) 

Figure B2 shows the three coefficients for an isotropic 
case as a function of the effective optical depth of the 
medium r e . As limiting values, we obtain 

a-jlTx" V ■ IBI7) 





and 



(tl 

1 (r L 



(r L ^ 


00) 


(r L ^ 


0) 


(tl- 


00) 


(r L ^ 


0) 


-> 00) 




-0) 





1 + x , (B18) 





(B19) 



